version 17.0
clear all
cd "MYPATH\descript_by_xs"
adopath + ../ado/
pause on
cap log close


log using analysis.log, replace

preliminaries
foreach PATH in RESULTS TEMP {
	cap mkdir "${`PATH'}\descript_by_xs"
    if "`PATH'" != "TEMP" cap mkdir "${`PATH'}\descript_by_xs\figures"
	if "`PATH'" != "TEMP" cap mkdir "${`PATH'}\descript_by_xs\fig_data"
}

graph set window fontface default
graph set ps fontface default
graph set window fontfacemono "Consolas"
graph set ps fontfacemono "Consolas"

program main
	prep_data
	grad_by_xs
	reg_on_xs
	export_figures
	cap log close
end


program prep_data
	foreach gr in all thou {

		if "`gr'" == "all" use "${DATA}\all_screens_universalnt_w_xs.dta", clear
		if "`gr'" == "thou" use "${DATA}\model_sample.dta", clear
		
		* age bins
		qui gen age_bin = . 
		forval i = 20(3)40 {
			replace age_bin = `i' + 1 if age >= `i' & age < (`i' + 3)
		}
		replace age_bin = 38 if age == 40
		replace age_bin = 41 if age >= 40
		replace age_bin = 19 if age < 20
		tab age_bin
		label def age_bin_lab 19 "< 20" 21 "20-22" 24 "23-25" 27 "26-28" 30 "29-31" 33 "32-34" 36 "35-37" ///
		  39 "38-40" 41 "> 40"
		label val age_bin age_bin_lab
		local decile = 10 
		gen inc_decile = .
		forval i = 0(10)100 {
			replace inc_decile = `i' + 5 if inc_rank >= `i' & inc_rank < (`i' + 10)
		}
		replace inc_decile = 95 if inc_decile >= 95
		
		gen age_bin_broad = .
		replace age_bin_broad = 1 if age < 25
		replace age_bin_broad = 2 if age >= 25 & age < 35
		replace age_bin_broad = 3 if age >= 35
		assert !mi(age_bin_broad)
		
		label def inc_quartilelab 1 "Income: q1" 2 "Income: q2" 3 "Income: q3" 4 "Income: q4"
		label def age_bin_broadlab 1 "Age < 25" 2 "Age between 25 and 35" 3 "Age > 35"
		label val inc_quartile inc_quartilelab 
		label val age_bin_broad age_bin_broadlab
		label var married "Married"
		label var mom_foreign "Foreign-born"
		label var dv_prev_kids "Any previous kids"
		label var prev_birth_issue "Prev birth issue"
		
		* KUB risk bins
		local bin_size = 50
		cap drop bin_number
		qui gen bin_number = .
		local count = 1
		forval bin_max = `bin_size'(`bin_size')20000 {
			local bin_inf = `bin_max' - `bin_size'
			qui replace bin_number = `count' if fetus_risk <= `bin_max' & fetus_risk > `bin_inf'
			local count = `count' + 1
		}

		foreach var in subt did_nipt did_invasive {
			replace `var' = `var' * 100
			qui sum `var'
			local mean_`var' = r(mean)
			qui areg `var', absorb(bin_number)
			predict double resid_`var', resid
			replace resid_`var' = resid_`var' + `mean_`var''
		}
		save "$TEMP\descript_by_xs\by_xs_`gr'_sample", replace
	}
end



program grad_by_xs
	foreach gr in all thou {
		use "$TEMP\descript_by_xs\by_xs_`gr'_sample", replace
		local share_vars subt did_nipt did_invasive
		*local mom_vars inc_rank age
		local mom_vars inc_rank age educ married mom_foreign dv_prev_kids prev_birth_issue age_bin_broad inc_quartile
		local share_vars_resids resid_subt resid_did_nipt resid_did_invasive	
		foreach mvar in `mom_vars' {
			local bins = "`mvar'"
			if "`mvar'" == "inc_rank" { 
				local bins = "inc_decile"
				local xsc = "0 100"
				local xlab = "0(10)100"
				local xtit = "Household income rank"
				local ylab = "0(2)10"
				local yscl = "0 10"
			}
			if "`mvar'" == "age" {
				local bins = "age_bin"
				local xsc = "18 41"
				local xlab = "19 21 24 27 30 33 36 39 41, valuelabel"
				local xtit = "Mother's age"
				local ylab = "0(5)20"
				local yscl = "0 20"
			}
			di "`mvar'"
			di "Number not missing `mvar'"
			count if !mi(`mvar')
			* Unconditional graphs
			preserve
			collapse (mean) `share_vars' `share_vars_resids', by(`bins')
			if "`mvar'" == "age" | "`mvar'" == "inc_rank" {
				twoway (connect subt `bins', m(o) mcolor(green) lcolor(green)) ///
				  (connect did_nipt `bins', m(t) mcolor(blue) lcolor(blue)) ///
				  (connect did_invasive `bins', m(d) mcolor(red) lcolor(red)), ///
				  name(grad_`mvar'_`gr') ///
				  legend(order(1 "Any post-NT testing" 2 "cfDNA" 3 "Invasive") ///
					pos(6) span row(1) col(3)) ///
				  xscale(r(`xsc')) xlab(`xlab') xtitle(`xtit', size(small)) ///
				  yscale(r(`yscl')) ytitle("" , size(small)) ylabel(`ylab', format(%3.0f))
			
				  
				* Conditional graphs
				twoway (connect resid_subt `bins', m(o) mcolor(green) lcolor(green)) ///
				  (connect resid_did_nipt `bins', m(t) mcolor(blue) lcolor(blue)) ///
				  (connect resid_did_invasive `bins', m(d) mcolor(red) lcolor(red)), ///
				  name(grad_cond_`mvar'_`gr') ///
				  legend(order(1 "Any post-NT testing" 2 "cfDNA" 3 "Invasive") ///
					pos(6) span row(1) col(3)) ///
				  xscale(r(`xsc')) xlab(`xlab') xtitle(`xtit', size(small)) ///
				  yscale(r(`yscl')) ytitle("" , size(small)) ylabel(`ylab', format(%3.0f))
			}
			local shares `share_vars' `share_vars_resids'
			foreach var of local shares {
				replace `var' = `var'/100
			}
			export delimited using "$RESULTS/descript_by_xs/fig_data/gradient_`mvar'_`gr'_new.csv", replace
			restore  
		}
		
		local mom_vars_discrete educ married mom_foreign dv_prev_kids prev_birth_issue age_bin_broad inc_quartile
		label_vals
		foreach mvar in `mom_vars_discrete' {
			di "`mvar'"
			di "Number not missing `mvar'"
			count if !mi(`mvar')
			* Unconditional graphs
			preserve
			collapse (mean) `share_vars' `share_vars_resids' (count) pregnancy, by(`mvar')
			rename pregnancy n_preg
			expand 3
			bysort `mvar': gen test_type = _n 
			gen share = subt if test_type == 1
			gen share_resid = resid_subt if test_type == 1
			replace share = did_nipt if test_type == 2
			replace share_resid = resid_did_nipt if test_type == 2
			replace share = did_invasive if test_type == 3
			replace share_resid = resid_did_invasive if test_type == 3
			label def type_lab 1 "Any post-NT testing" 2 "cfDNA" 3 "Invasive"
			label val test_type type_lab 
			drop `share_vars' `share_vars_resids'
			graph bar share, over(`mvar') over(test_type) asyvars name(grad_`mvar'_`gr') ///
			  bar(1, color(green)) bar(2, color(blue)) bar(3, color(red)) ///
			  blabel(total, format(%3.2fc)) ///
			  legend(pos(6) span row(1)) ytitle("")
			graph bar share_resid, over(`mvar') over(test_type) asyvars name(grad_cond_`mvar'_`gr') ///
			  bar(1, color(green)) bar(2, color(blue)) bar(3, color(red)) ///
			  blabel(total, format(%3.2fc)) ///
			  legend(pos(6) span row(1)) ytitle("")
			export delimited using "$RESULTS/descript_by_xs/fig_data/gradient_`mvar'_`gr'.csv", replace
			restore 
		}
	}
end

program label_vals
	label def marrylab 0 "Not married" 1 "Married"
	label val married marrylab
	label def foreignlab1 0 "Swedish" 1 "Foreign-born"
	label val mom_foreign foreignlab1
	label def dv_prev_kidslab 0 "No prev kids" 1 "Prev kids"
	label val dv_prev_kids dv_prev_kidslab
	label def prev_birth_issuelab 0 "No prev birth issue" 1 "Prev birth issue"
	label val prev_birth_issue prev_birth_issuelab 
end

program reg_on_xs
	foreach gr in all thou {
		use "$TEMP\descript_by_xs\by_xs_`gr'_sample", clear
		drop age_bin
		qui gen age_bin = 1 if (age < 25)
		qui replace age_bin = 2 if (age >= 25 & age < 35)
		qui replace age_bin = 3 if (age >= 35 & !mi(age))
		label var age_bin "Age" 
		label def agelab 1 "< 25" 2 "25-30" 3 "35+"
		label val age_bin agelab
		replace edu = 4 if (mi(educ))
		label def labeduc 1 "No college" 2 "Some college" 3 "Full college" 4 "Missing"
		label val educ labeduc
		replace inc_quartile = 5 if (mi(inc_quartile))
		label def incquartilelab 1 "Income: q1" 2 "Income: q2" 3 "Income: q3" /// 
		4 "Income: q4" 5 "Missing"
		label val inc_quartile incquartilelab 

		local share_vars_regs subt resid_subt did_nipt resid_did_nipt did_invasive resid_did_invasive
		foreach var in `share_vars_regs' {
			qui eststo `var'_reg_xs: reg `var' i.age_bin i.edu ///
			  i.inc_quartile ///
			  married mom_foreign dv_prev_kids ///
			  prev_birth_issue, robust   
		}
		
		foreach svar in subt did_nipt did_invasive {
			foreach var in age_bin edu inc_quartile married mom_foreign dv_prev_kids prev_birth_issue age_bin_broad {
				local xvar = "`var'"
				if inlist("`var'", "age_bin", "edu", "inc_quartile") local xvar = "i.`var'"
				if "`svar'" == "subt" local s = "subt"
				if "`svar'" == "did_invasive" local s = "inv"
				if "`svar'" == "did_nipt" local s = "nipt"
				qui eststo `s'_reg_`var': reg `svar' `xvar', robust   
			}
		}		
		
		qui cd "$RESULTS\descript_by_xs\"
		esttab subt_reg_xs resid_subt_reg_xs did_nipt_reg_xs resid_did_nipt_reg_xs ///
		  did_invasive_reg_xs resid_did_invasive_reg_xs ///
		  using ".\testing_reg_xs_`gr'.tex", ///
		  stats(N) se label replace booktabs ///
		  mtitles("Unconditional" "Conditional" "Unconditional" "Conditional" "Unconditional" "Conditional")
		  
		esttab subt_reg_xs resid_subt_reg_xs did_nipt_reg_xs resid_did_nipt_reg_xs ///
		  did_invasive_reg_xs resid_did_invasive_reg_xs ///
		  using ".\testing_reg_xs_`gr'.csv", ///
		  stats(N) se label replace ///
		  mtitles("Unconditional" "Conditional" "Unconditional" "Conditional" "Unconditional" "Conditional")
	 
		esttab subt_reg_age_bin subt_reg_edu subt_reg_inc_quartile ///
		  subt_reg_married subt_reg_mom_foreign ///
		  subt_reg_dv_prev_kids subt_reg_prev_birth_issue ///
		  using ".\subt_reg_indiv_xs_`gr'.tex", ///
		  stats(N) se label replace booktabs 
		
		esttab inv_reg_age_bin inv_reg_edu inv_reg_inc_quartile ///
		  inv_reg_married inv_reg_mom_foreign ///
		  inv_reg_dv_prev_kids inv_reg_prev_birth_issue ///
		  using ".\inv_reg_indiv_xs_`gr'.tex", ///
		  stats(N) se label replace booktabs 
		
		esttab nipt_reg_age_bin nipt_reg_edu nipt_reg_inc_quartile ///
		  nipt_reg_married nipt_reg_mom_foreign ///
		  nipt_reg_dv_prev_kids nipt_reg_prev_birth_issue ///
		  using ".\nipt_reg_indiv_xs_`gr'.tex", ///
		  stats(N) se label replace booktabs 
	}
end

program export_figures
	foreach gr in all thou {
		foreach mvar in inc_rank age educ married mom_foreign dv_prev_kids prev_birth_issue age_bin_broad inc_quartile {
			di "`mvar', `gr'"
			graph export "${RESULTS}\descript_by_xs\figures/grad_`mvar'_`gr'.pdf", as(pdf) replace ///
				name(grad_`mvar'_`gr')
			graph export "${RESULTS}\descript_by_xs\figures/grad_cond_`mvar'_`gr'.pdf", as(pdf) replace ///
				name(grad_cond_`mvar'_`gr')
		}
	}
end



* Execute 
main
